%%
iter_PRIVATE=0;
iter_tol_PRIVATE=1000;
DPRIVATE=10;
tol_PRIVATE=1e-7;
outfreq_PRIVATE=250;

%For the Next Loop

bp_cc=bp;%zeros(NBG,NB,NSS,NBG);

iter_public=0;
outfreq_public=10;
outfreq_public_private=50;
dist_public=10;
maxit_public=500;
tol_public=1e-5;
uptd=0.2;
%      evdef_big_2d=evdef_big_2d%squeeze(evdef_big(1,:,:));
      bp_def_2d=bp_bad; %squeeze(bp_def(1,:,:));

%%    A Peusdo Private Equilibruim

while dist_public>tol_public && iter_public<maxit_public
     


PreQ=(1/R)*Repay.*(delta*ones(NBG,NB,NSS)+(1-delta)*sum(Qgov.*G_big,4));


             parfor j=1:NSS
                        Ptemp=Prob(j,:)';
                        bdef_badS=squeeze(bp_def_2d(:,j)); % THIS IS MY BORROWING CHOICE IN BAD STANDING
                        evdef_temp_badS=VbadS*Ptemp; % THIS IS THE BAD STANDING VALUE
                        vdef_badS_interp=griddedInterpolant(bgrid1d,evdef_temp_badS,'linear','linear');
                        evdef_badS(:,j)=vdef_badS_interp(bdef_badS);     % EXPECTED VALUE IF STAY IN BAD STANDING  
                        reentry_temp=squeeze(V(zeroBG_ind,:,:))*Ptemp; % THIS IS THE REENTRY VALUE
                        vreentry_interp=griddedInterpolant(bgrid1d,reentry_temp,'linear','linear');
                        evdef(:,j)=vreentry_interp(bdef_badS); % EXPECTED VALUE IF REENTERS
                        for inext=1:NBG
                          evcc_temp=squeeze(V(inext,:,:))*Ptemp;
                          vcc_interp=griddedInterpolant(bgrid1d,evcc_temp,'linear','linear');
                          bc0=squeeze(bp(:,:,j,inext));
                          bc=bc0(:);
                          evcBigR(:,j,inext)=vcc_interp(bc);
                       end
                end 
                evcBig=reshape(evcBigR(:),NBG,NB,NSS,NBG);

     
            
      parfor j=1:NSS      
          Ptemp=Prob(j,:)';
        evdef_temp=squeeze(V(zeroBG_ind,:,:))*Ptemp; % VALUE OF REENTRY
        evdef_temp_badS=VbadS*Ptemp; % THIS IS THE BAD STANDING VALUE
        vdef_interp=griddedInterpolant(bgrid1d,evdef_temp,'linear','linear');
        vdef_badS_interp=griddedInterpolant(bgrid1d,evdef_temp_badS,'linear','linear');
        evdef_temp_large(:,j)=vdef_interp(squeeze(bgrid_large(:))); %REENTRY
        evdef_temp_large_badS(:,j)=vdef_badS_interp(squeeze(bgrid_large(:))); %BAD STANDING
        for inext=1:NBG
          evcc_temp=squeeze(V(inext,:,:))*Ptemp;
          eQ_temp=squeeze(PreQ(inext,:,:))*Ptemp;
          vcc_interp=griddedInterpolant(bgrid1d,evcc_temp,'linear','linear');
          eq_interp=griddedInterpolant(bgrid1d,eQ_temp,'linear','linear');
          evcc_temp_large(inext,:,j)=vcc_interp(bgrid_large(:));
          eQ_temp_large(inext,:,j)=eq_interp(bgrid_large(:));
        end
      end   
        eQ_temp_large(eQ_temp_large<0)=1e-10;
        eQ_temp_large(eQ_temp_large>delta/(rbar+delta))=delta/(rbar+delta);
        oldQ=Qgov;
        oldevc=evc_big;
        old_bp=bp_cc;

        evdef_final=theta*evdef_temp_large+(1-theta)*evdef_temp_large_badS;
      old_evdef_big_2d=evdef_big_2d; 
      old_bp_def_2d=bp_def_2d;
      
%  The multiple loops

   parfor j=1:NSS %Today's Exogenous Shock
        %Endownment
        ytemp=yT(j);
        yntemp=yN(j);
        qpv=q(j);
        pvd_temp=PVD(j);
        kappa_temp=KAPPAS(j);

           for t=1:NB %Today's private debt index                 
               %Compute Vd
                        btemp=bgrid1d(t);                      
                        current_c_def=1e-8;
                        current_evdef=old_evdef_big_2d(t,j)
                        current_max_def=((omega*current_c_def^(-ita)+(1-omega)*yntemp^(-ita))^((sigma-1)/ita)-1)/(1-sigma)-NUS(j)+beta*current_evdef;
                        current_bp_def=old_bp_def_2d(t,j);
                        ups_def_temp=0;
                        for x=1:NB2
                        bgdtc_def=0;
                        bnext_temp_def=bgrid_large(x);                        
                        evdef=evdef_final(x,j);
                        c_def_temp=max(ytemp+bnext_temp_def*qpv-(1-pvd_temp)*btemp,1e-8); 
                        price_def_temp=real((1-omega)/omega*(c_def_temp/yntemp).^(1+ita));
                        ups_def_temp=((omega*c_def_temp^(-ita)+(1-omega)*yntemp^(-ita))^((sigma-1)/ita)-1)/(1-sigma)-NUS(j)+beta*evdef;
                        if bnext_temp_def-(1/qpv)*kappa_temp*(price_def_temp*yntemp + ytemp) >1e-8 && c_def_temp>1e-8
                             bgdtc_def=1; %-250
                        end
                            if ups_def_temp>current_max_def && bgdtc_def==0
                                current_max_def=ups_def_temp;
                                current_c_def=c_def_temp;
                                current_bp_def=bnext_temp_def;
                                current_evdef=evdef;
                            end                       
                        end
                        Ups_d_2d(t,j)=current_max_def;
                        c_def_2d(t,j)=current_c_def;
                        bp_def_2d(t,j)=current_bp_def;
                        evdef_big_2d(t,j)=current_evdef;
                        
                for i=1:NBG %Today's public Debt shock
                     BGtemp=BGgrid(i);    
                        for inext=1:NBG %This is the Grid Search                            
                            BGnext_temp=BGgrid(inext);
                            current_c_cc=1e-8;
                            current_bp_cc=old_bp(i,t,j,inext);
                            qtemp=oldQ(i,t,j,inext);
                            current_evc=oldevc(i,t,j,inext);
                            current_max_cc=((omega*current_c_cc^(-ita)+(1-omega)*yntemp^(-ita))^((sigma-1)/ita)-1)/(1-sigma)+beta*current_evc;
                            current_q=qtemp;
                            current_T=current_q*(BGnext_temp-(1-delta)*BGtemp)-delta*BGtemp;
                            ups_cc_temp=0;
                          for xx=1:NB2  
                              bgdtc_cc=0;
                            bnext_temp=bgrid_large(xx);
                            qtemp=eQ_temp_large(inext,xx,j);
                            evc=evcc_temp_large(inext,xx,j);
                            c_c_big_temp=max(ytemp+bnext_temp*qpv-(1-pvd_temp)*btemp+qtemp*(BGnext_temp-(1-delta)*BGtemp)-delta*BGtemp,1e-8);                            
                            price_cc_temp=real((1-omega)/omega*(c_c_big_temp/yntemp)^(1+ita));
                            ups_cc_temp=((omega*c_c_big_temp^(-ita)+(1-omega)*yntemp^(-ita))^((sigma-1)/ita)-1)/(1-sigma)+beta*evc;
                            if bnext_temp-(1/qpv)*kappa_temp*(price_cc_temp*yntemp + ytemp)>1e-8 && c_c_big_temp>1e-8
                                bgdtc_cc=1; %-250                           
                            end
                          
                            if ups_cc_temp>current_max_cc && bgdtc_cc==0
                                current_max_cc=ups_cc_temp;
                                current_c_cc=c_c_big_temp;
                                current_q=qtemp;
                                current_T=qtemp*(BGnext_temp-(1-delta)*BGtemp)-delta*BGtemp;
                                %current_bpi_cc=xx;
                                current_bp_cc=bnext_temp;
                                current_evc=evc;
                            end   
                          end
                          if current_c_cc>1e-8
                            Ups_big(i,t,j,inext)=current_max_cc;
                            Qgov(i,t,j,inext)= current_q; 
                            c_c_big(i,t,j,inext)=current_c_cc;
                            Transfer(i,t,j,inext)=current_T;
                            bp_cc(i,t,j,inext)=current_bp_cc;
                            evc_big(i,t,j,inext)=current_evc;
                          else % AS IN THE BASELINE WHEN THE CREDIT CONSTRAINT CANNOT HOLD WITH POSITIVE CONSUMPTION WE ALLOW FOR VIOLATION IT
                           c_c_big(i,t,j,inext)=current_c_cc;
                           current_bp_cc=(1/qpv)*(current_c_cc-(ytemp-(1-pvd_temp)*btemp+qtemp*(BGnext_temp-(1-delta)*BGtemp)-delta*BGtemp));
                            bp_cc(i,t,j,inext)=current_bp_cc; 
                            Qterp=griddedInterpolant(squeeze(bgrid_large(:)),squeeze(eQ_temp_large(inext,:,j)),'linear','linear');
                            Vterp=griddedInterpolant(squeeze(bgrid_large(:)),squeeze(evcc_temp_large(inext,:,j)),'linear','linear');
                            current_q=Qterp(current_bp_cc);
                            current_evc=Vterp(current_bp_cc);                           
                            Qgov(i,t,j,inext)= current_q; 
                            Transfer(i,t,j,inext)=current_q*(BGnext_temp-(1-delta)*BGtemp)-delta*BGtemp;
                            evc_big(i,t,j,inext)=current_evc;
                            Ups_big(i,t,j,inext)=((omega*current_c_cc^(-ita)+(1-omega)*yntemp^(-ita))^((sigma-1)/ita)-1)/(1-sigma)+beta*current_evc;                           

                          end   
                        end 

                 end
           end
   end
 % RESIZE DEFAULT
  ups_def_temp=((omega*c_def_2d.^(-ita)+(1-omega)*yN2D.^(-ita)).^((sigma-1)/ita)-1*ones(NB,NSS))/(1-sigma)-nus2D+beta*(evdef_big_2d);                  
  Ups_d=reshape(repmat(ups_def_temp(:),1,NBG)',NBG,NB,NSS);

   %Repayment

               
                 %Continuation
                 Ups_big=((omega*c_c_big.^(-ita)+(1-omega)*yn_big.^(-ita)).^((sigma-1)/ita)-1*ones(NBG,NB,NSS,NBG))/(1-sigma)+beta*evc_big;
                 Ups_big(c_c_big==1e-8)=-2e6;             
                 %Update Probabilities
                 Upsmax_small=max(Ups_big,[],4);
                 Upsmax = repmat(Upsmax_small,1,1,1,NBG); 
                 %WE ASSUME CORRELATED CHOICES 
                sumUR_small=(sum(exp((Ups_big-Upsmax)/(v*p)),4)).^(p);
                %CONDITIONAL PROBABILITES UNDER REPAYMENT
                G_big=exp((Ups_big-Upsmax)/(v*p))./repmat(sumUR_small.^(1/p),1,1,1,NBG);
                % TOTAL MAX
                UpsmaxTot=max(Upsmax_small,Ups_d);
                %COMPUTE V
                IVALUE0=exp((Ups_d-UpsmaxTot)/v)+exp((Upsmax_small-UpsmaxTot)/v).*sumUR_small;
                V=UpsmaxTot+v*log(IVALUE0);
                VbadS=ups_def_temp+mean_logi_shocks*ones(NB,NSS); % THE VALU IN BAD STANDING IS UPS_DEF + THE AVERAGE SHOCK BECAUSE IID OVER TIME
                % PROBABILITY OF DEFAULT and Repayment
                Def_func=exp((Ups_d-UpsmaxTot)/v)./IVALUE0;
                Repay=(exp((Upsmax_small-UpsmaxTot)/v).*sumUR_small)./IVALUE0;
                % UNCONDITIONAL PROBABILITIES
                BigRepay=repmat(Repay,1,1,1,NBG);


        % Update Prices
          
     
        dist_public=squeeze(squeeze(max(max(max(abs(V_upd-V))))));
        V_upd=V;
        iter_public=iter_public+1;
        if mod(iter_public, outfreq_public) == 0 || iter_public==1
            disp('Public Outer      ');
            fprintf('%d          %1.7f \n',iter_public,dist_public);
        end
        if dist_public<10*tol_public
            outfreq_public_private=1;
            outfreq_public=1;
        end
        tV=isnan(V);        
        tQ=isnan(Q);        
        
        if max(tV(:))>0 ||max(tQ(:))>0 %||max(max(max(tEmupx)))>0||max(max(max(mUps)))>0
            stop
        end
       

end



